{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian 中的 PUHF/PMP2 结果的重新实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2019-08-31，最后修改：2019-09-01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这一份笔记中，我们将使用 PySCF 的功能与 NumPy 重复 Gaussian 中计算的 PUHF 与 PMP2 能量结果；并对 PUHF 与 PMP2 的推导作简单的说明。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf, mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考结果与体系定义"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian 结果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 Gaussian 中，我们使用以下输入卡可以得到 PUHF/PMP2 能量：\n",
    "\n",
    "```\n",
    "#p UMP2(Full)/6-31G nosymm\n",
    "\n",
    "H2O\n",
    "\n",
    "3 4\n",
    "O  0.  0.  0.\n",
    "H  1.  0.  0.\n",
    "H  0.  1.  0.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于上述分子，其中一些重要的输出结果是\n",
    "\n",
    "1. $E_\\mathrm{UHF}$：-73.0451423839\n",
    "\n",
    "2. $E_\\mathrm{UMP2, corr}$: -0.02646719276\n",
    "\n",
    "3. $E_\\mathrm{UMP2}$: -73.071609576661\n",
    "\n",
    "4. $\\langle S_z \\rangle$: 1.5\n",
    "\n",
    "5. $\\langle S^{2(0)} \\rangle$: 3.7531\n",
    "\n",
    "6. $\\langle S^{2(0)} \\rangle + \\langle S^{2(1)} \\rangle$: 3.7504\n",
    "\n",
    "7. $E_\\mathrm{PUHF}$：-73.046146318\n",
    "\n",
    "8. $E_\\mathrm{PMP2}$: -73.072180589\n",
    "\n",
    "输出文件参见 {download}`assets/PUHF_and_PMP2.out`；其中，有效的数据可以通过下述的代码获得："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "line 336:  SCF Done:  E(UHF) =  -73.0451423839     A.U. after   11 cycles\n",
      "line 338:  <Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 1.5000 <S**2>= 3.7531 S= 1.5008\n",
      "line 370:  (S**2,0)=  0.37531D+01           (S**2,1)=  0.37504D+01\n",
      "line 371:  E(PUHF)=      -0.73046146318D+02        E(PMP2)=      -0.73072180589D+02\n",
      "line 373:  E2 =    -0.2646719276D-01 EUMP2 =    -0.73071609576661D+02\n"
     ]
    }
   ],
   "source": [
    "with open(\"assets/PUHF_and_PMP2.out\", \"r\") as output:\n",
    "    output_lines = output.read().split(\"\\n\")\n",
    "\n",
    "for line_num, line_text in enumerate(output_lines):\n",
    "    if any([keyword in line_text for keyword in\n",
    "            [\"SCF Done\", \"EUMP2\", \"<S**2>\", \"(S**2,1)\", \"E(PMP2)\"]]) \\\n",
    "        and \"Initial guess\" not in line_text:\n",
    "        print(\"line {:03d}: {}\".format(line_num, line_text))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们的目标就是近乎准确无误地重复上述八个结果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PySCF 体系定义"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了获得与 Gaussian 相同的结果，我们需要定义相同的分子与电荷、多重度环境："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7fa4e5342160>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "O 0. 0. 0.\n",
    "H 1. 0. 0.\n",
    "H 0. 1. 0.\n",
    "\"\"\"\n",
    "mol.charge = 3\n",
    "mol.spin = 3\n",
    "mol.basis = \"6-31G\"\n",
    "mol.build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过 PySCF 计算 UHF 能量："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged SCF energy = -73.0451423536459  <S^2> = 3.7530824  2S+1 = 4.0015409\n"
     ]
    }
   ],
   "source": [
    "scf_eng = scf.UHF(mol)\n",
    "scf_eng.conv_tol = 1e-10\n",
    "scf_eng.run();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述结果应当能与 $E_\\mathrm{UHF}$ 和 $\\langle S^{2(0)} \\rangle$ 对应。$\\langle S_z \\rangle = 1.5$ 几乎是显然的。不过，我们仍然不了解 $\\langle S^{2(0)} \\rangle$ 是如何生成的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过 PySCF 计算 UMP2 能量："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(UMP2) = -73.0716095455079  E_corr = -0.0264671918619837\n"
     ]
    }
   ],
   "source": [
    "mp2_eng = mp.UMP2(scf_eng)\n",
    "mp2_eng.run();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述结果应当能与 $E_\\mathrm{UMP2, corr}$ 和 $E_\\mathrm{UMP2}$ 对应。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此，当前的问题将是回答：如何重复\n",
    "\n",
    "1. $\\langle S^{2(0)} \\rangle$: 3.7531\n",
    "\n",
    "2. $\\langle S^{2(0)} \\rangle + \\langle S^{2(1)} \\rangle$: 3.7504\n",
    "\n",
    "3. $E_\\mathrm{PUHF}$：-73.046146318\n",
    "\n",
    "4. $E_\\mathrm{PMP2}$: -73.072180589"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 部分变量定义"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先，我们遵从大多数量化文章中的记号\n",
    "\n",
    "- $i, j$ 代表占据分子轨道\n",
    "\n",
    "- $a, b$ 代表非占分子轨道\n",
    "\n",
    "- $p, q$ 代表任意分子轨道\n",
    "\n",
    "- $\\alpha, \\beta$ 代表任意原子轨道"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><b>Table 1. 分子相关变量</b></center>\n",
    "\n",
    "| 变量名 | 元素记号 | 意义与注解 | 标量或区间 |\n",
    "|-|-|-|-|\n",
    "| `nocc_a` | $n_\\mathrm{occ}^\\alpha$ | $\\alpha$ 自旋电子数 | $5$ |\n",
    "| `nocc_b` | $n_\\mathrm{occ}^\\beta$ | $\\beta$ 自旋电子数 | $2$ |\n",
    "| `N` | $N$ | 总电子数 | $7$ |\n",
    "| `nmo` | $n_\\mathrm{MO}$ | 分子轨道数 | $13$ |\n",
    "| `nao` | $n_\\mathrm{AO}$ | 原子轨道数 | $13$ |\n",
    "| `S` | $S_{\\mu \\nu}$ | 原子轨道重叠积分 | |\n",
    "| `so_a` | | $\\alpha$ 占据轨道分割 | $[0, 5)$ |\n",
    "| `so_b` | | $\\beta$ 占据轨道分割 | $[0, 2)$ |\n",
    "| `sv_a` | | $\\alpha$ 非占轨道分割 | $[5, 13)$ |\n",
    "| `sv_b` | | $\\beta$ 非占轨道分割 | $[2, 13)$ |\n",
    "| `Sx` | $S_x$ | $x$ 分量自旋 | $0$ |\n",
    "| `Sy` | $S_y$ | $y$ 分量自旋 | $0$ |\n",
    "| `Sz` | $S_z$ | $z$ 分量自旋 | $3/2$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><b>Table 2. UHF 计算相关变量</b></center>\n",
    "\n",
    "| 变量名 | 元素记号 | 意义与注解 |\n",
    "|-|-|-|\n",
    "| `C_a` | $C_{\\mu p}$ | $\\alpha$ 系数矩阵 |\n",
    "| `C_b` | $C_{\\mu \\bar p}$ | $\\beta$ 系数矩阵 |\n",
    "| `e_a` | $e_{p}$ | $\\alpha$ 轨道能 |\n",
    "| `e_b` | $e_{\\bar p}$ | $\\alpha$ 轨道能 |\n",
    "| `eo_a` | $e_{i}$ | $\\beta$ 占据轨道能 |\n",
    "| `eo_b` | $e_{\\bar i}$ | $\\alpha$ 占据轨道能 |\n",
    "| `ev_a` | $e_{a}$ | $\\alpha$ 非占轨道能 |\n",
    "| `ev_b` | $e_{\\bar a}$ | $\\beta$ 非占轨道能 |\n",
    "| `D2_aa` | $D_{ij}^{ab}$ | $\\alpha, \\alpha$ 轨道能差 |\n",
    "| `D2_ab` | $D_{i \\bar j}^{a \\bar b}$ | $\\alpha, \\beta$ 轨道能差 |\n",
    "| `D2_bb` | $D_{\\bar i \\bar j}^{\\bar a \\bar b}$ | $\\beta, \\beta$ 轨道能差 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><b>Table 3. UMP2 计算相关变量</b></center>\n",
    "\n",
    "| 变量名 | 元素记号 | 意义与注解 |\n",
    "|-|-|-|\n",
    "| `t2_aa` | $t_{ij}^{ab}$ | $\\alpha, \\alpha$ MP2 激发系数 |\n",
    "| `t2_ab` | $t_{i \\bar j}^{a \\bar b}$ | $\\alpha, \\beta$ MP2 激发系数 |\n",
    "| `t2_bb` | $t_{\\bar i \\bar j}^{\\bar a \\bar b}$ | $\\beta, \\beta$ MP2 激发系数 |\n",
    "| `D2_aa` | $D_{ij}^{ab}$ | $\\alpha, \\alpha$ MP2 激发系数分母 |\n",
    "| `D2_ab` | $D_{i \\bar j}^{a \\bar b}$ | $\\alpha, \\beta$ MP2 激发系数分母 |\n",
    "| `D2_bb` | $D_{\\bar i \\bar j}^{\\bar a \\bar b}$ | $\\beta, \\beta$ MP2 激发系数分母 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述需要补充说明的公式有：\n",
    "\n",
    "$$\n",
    "S_z = \\frac{1}{2} (n_\\mathrm{occ}^\\alpha - n_\\mathrm{occ}^\\beta)\n",
    "$$\n",
    "\n",
    "$$\n",
    "D_{i \\bar j}^{a \\bar b} = e_i + e_{\\bar j} - e_a - e_{\\bar b}\n",
    "$$\n",
    "\n",
    "对于 MP2 激发系数分母，另外两种自旋情况的 $D_{ij}^{ab}$ 与 $D_{\\bar i \\bar j}^{\\bar a \\bar b}$ 也可以类似地生成。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# === Molecular\n",
    "# --- Definition\n",
    "nocc_a, nocc_b = mol.nelec\n",
    "N = nocc_a + nocc_b\n",
    "nmo = nao = mol.nao\n",
    "S = mol.intor(\"int1e_ovlp\")\n",
    "# --- Derivative\n",
    "so_a, so_b = slice(0, nocc_a), slice(0, nocc_b)\n",
    "sv_a, sv_b = slice(nocc_a, nmo), slice(nocc_b, nmo)\n",
    "Sx, Sy, Sz = 0, 0, 0.5 * (nocc_a - nocc_b)\n",
    "\n",
    "# === UHF Calculation\n",
    "# --- Definition\n",
    "C_a, C_b = scf_eng.mo_coeff\n",
    "e_a, e_b = scf_eng.mo_energy\n",
    "# --- Derivative\n",
    "eo_a, eo_b = e_a[so_a], e_b[so_b]\n",
    "ev_a, ev_b = e_a[sv_a], e_b[sv_b]\n",
    "D2_aa = eo_a[:, None, None, None] + eo_a[None, :, None, None] - ev_a[None, None, :, None] - ev_a[None, None, None, :]\n",
    "D2_ab = eo_a[:, None, None, None] + eo_b[None, :, None, None] - ev_a[None, None, :, None] - ev_b[None, None, None, :]\n",
    "D2_bb = eo_b[:, None, None, None] + eo_b[None, :, None, None] - ev_b[None, None, :, None] - ev_b[None, None, None, :]\n",
    "\n",
    "# === MP2 Calculation\n",
    "t2_aa, t2_ab, t2_bb = mp2_eng.t2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作为对四脚标张量性质的验证，我们计算 MP2 相关能 $E_\\mathrm{MP2, corr}$ 如下：\n",
    "\n",
    "$$\n",
    "E_\\mathrm{MP2, corr} =\n",
    "\\frac{1}{4} \\sum_{ijab} (t_{ij}^{ab})^2 D_{ij}^{ab} +\n",
    "\\frac{1}{4} \\sum_{\\bar i \\bar j \\bar a \\bar b} (t_{\\bar i \\bar j}^{\\bar a \\bar b})^2 D_{i \\bar j}^{a \\bar b} +\n",
    "\\sum_{i \\bar j a \\bar b} (t_{i \\bar j}^{a\\bar b})^2 D_{i \\bar j}^{a \\bar b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.02646719186198365"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(+ 0.25 * (t2_aa**2 * D2_aa).sum()\n",
    " + 0.25 * (t2_bb**2 * D2_bb).sum()\n",
    " + (t2_ab**2 * D2_ab).sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PySCF 所给出的 $E_\\mathrm{MP2, corr}$ 可以给出相同的结果："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.02646719186198366"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp2_eng.e_corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\langle S^2 \\rangle$ 相关计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 分子轨道基组重叠矩阵 `S_pq` $S_{p \\bar q}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "S_{p \\bar q} = \\sum_{\\mu \\nu} C_{\\mu p} S_{\\mu \\nu} C_{\\nu \\bar q}\n",
    "$$\n",
    "\n",
    "若用量子力学记号，上述矩阵元素可能表示为\n",
    "\n",
    "$$\n",
    "S_{p \\bar q} = \\int \\phi_p (\\boldsymbol{r}) \\phi_{\\bar q} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
    "$$\n",
    "\n",
    "注意上述的积分是空间坐标的积分，不包含自旋部分的积分。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(13, 13)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_pq = C_a.T @ S @ C_b\n",
    "S_pq.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们以后还会使用上述矩阵的占据-占据部分 `S_ij` $S_{i \\bar j}$、占据-非占部分 `S_ia` $S_{i \\bar a}$ 与非占-占据部分 `S_ai` $S_{a \\bar i} = S_{\\bar i a}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(5, 2), (5, 11), (8, 2)]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S_ij, S_ia, S_ai = S_pq[so_a, so_b], S_pq[so_a, sv_b], S_pq[sv_a, so_b]\n",
    "[S_ij.shape, S_ia.shape, S_ai.shape]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `S2_0` $\\langle S^{2(0)} \\rangle$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\langle S^{2(0)} \\rangle$ 在程序中通常写为 `<S^2>` 或 `<S**2>`。在 Gaussian 计算 PUHF 处，还写为 `(S**2,0)`。这意味着是 UHF 波函数的 $\\langle S^2 \\rangle_\\mathrm{UHF}$。相对地，UMP2 波函数给出的对 $\\langle S^2 \\rangle$ 的矫正将记作 $\\langle S^{2(1)} \\rangle$。\n",
    "\n",
    "参考 Chen and Schlegel [^Chen-Schlegel.JCP.1994.101] Table 1, $0 \\rightarrow 0$ 或等价地，Szabo and Ostlund [^Szabo-Ostlund.Dover.1996] eq (2.271)\n",
    "\n",
    "$$\n",
    "\\langle S^{2(0)} \\rangle = \\langle \\Psi_0 | \\hat S^2 | \\Psi_0 \\rangle = S_z (S_z + 1) + n_\\mathrm{occ}^\\beta - \\sum_{i \\bar j} (S_{i \\bar j})^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.7530823820890378"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S2_0 = Sz * (Sz + 1) + nocc_b - (S_ij**2).sum()\n",
    "S2_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian 的参考值是 3.7531。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了以后的记号便利，我们在这里定义 `L`\n",
    "\n",
    "$$\n",
    "L = \\sum_{i \\bar j} (S_{i \\bar j})^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = (S_ij**2).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `S2_1` $\\langle S^{2(1)} \\rangle$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\langle S^{2(1)} \\rangle &= 2 \\langle \\Psi_0 | \\hat S^2 | \\Psi^{(1)} \\rangle = 2 \\sum_{i \\bar j a \\bar b} t_{i \\bar j}^{a \\bar b} \\langle \\Psi_0 | \\hat S^2 | \\Psi_{i \\bar j}^{a \\bar b} \\rangle \\\\\n",
    "&= - 2 \\sum_{i \\bar j a \\bar b} t_{i \\bar j}^{a \\bar b} \\langle i | \\bar b \\rangle \\langle a | \\bar j \\rangle = - 2 \\sum_{i \\bar j a \\bar b} t_{i \\bar j}^{a \\bar b} S_{i \\bar b} S_{a \\bar j}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上式的第一个等号是 Chen and Schlegel [^Chen-Schlegel.JCP.1994.101] eq (5) 所给出的；而第三个等号是 Table 1 $0 \\rightarrow \\alpha \\beta (i, a: \\alpha; j, b: \\beta)$ 给出的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上式的推导中有一处关于 $| \\Psi^{(1)} \\rangle$ 的展开的推导省略。我们知道\n",
    "\n",
    "$$\n",
    "| \\Psi^{(1)} \\rangle = \\hat T_2 | \\Psi_0 \\rangle\n",
    "= \\frac{1}{4} \\sum_{ijab} t_{ij}^{ab} | \\Psi_{ij}^{ab} \\rangle + \\frac{1}{4} \\sum_{\\bar i \\bar j \\bar a \\bar b} t_{\\bar i \\bar j}^{\\bar a \\bar b} | \\Psi_{\\bar i \\bar j}^{\\bar a \\bar b} \\rangle + \\sum_{i \\bar j a \\bar b} t_{i \\bar j}^{a \\bar b} | \\Psi_{i \\bar j}^{a \\bar b} \\rangle\n",
    "$$\n",
    "\n",
    "但由于利用到 $\\langle 0 | \\hat S^2 | \\Psi_{ij}^{ab} \\rangle = \\langle 0 | \\hat S^2 | \\Psi_{\\bar i \\bar j}^{\\bar a \\bar b} \\rangle = 0$，因此在第二个等号时只从三个 $| \\Psi^{(1)} \\rangle$ 中留下了一项。关于 $\\hat S^2$ 作用在 UHF 波函数与轨道下的性质，可以参考 Schlegel [^Schlegel-Schlegel.JCP.1986.84] eq (5) 的说明。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.002658400959213991"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S2_1 = - 2 * (t2_ab * S_ia[:, None, None, :] * S_ai.T[None, :, :, None]).sum()\n",
    "S2_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此，UMP2 矫正过的 $\\langle S^2 \\rangle_\\mathrm{UMP2} = \\langle S^{2(0)} \\rangle + \\langle S^{2(1)} \\rangle$ 的结果是"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.7504239811298237"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S2_0 + S2_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian 的参考值是 3.7504。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `S4SD` $\\texttt{S4SD}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`S4SD` 的表达式较为复杂，我们也直接使用 $\\texttt{S4SD}$ 而不用其他记号表示该项：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\texttt{S4SD} = (n_\\mathrm{occ}^\\alpha - L) (n_\\mathrm{occ}^\\beta - L) + 2 L - 2 \\sum_{i \\bar j k \\bar l} S_{i \\bar j} S_{\\bar j k} S_{k \\bar l} S_{\\bar l i} + \\langle S^{2(0)} \\rangle^2\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.101029785519533"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S4SD = (nocc_a - L) * (nocc_b - L) + 2 * L - 2 * (S_ij @ S_ij.T @ S_ij @ S_ij.T).trace() + S2_0**2\n",
    "S4SD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该表达式的来源可能是 Amos and Hall [^Amos-Hall.PRSLA.1961.263]。该文的 eq (7·02) 下方公式中，有通过稍高阶的投影而获得的 $\\langle S^2 \\rangle$ 的计算方式\n",
    "\n",
    "$$\n",
    "\\langle S^2 \\rangle \\simeq \\langle S^{2(0)} \\rangle + \\frac{\\texttt{S4SD} - \\langle S^{2(0)} \\rangle^2}{\\langle S^{2(0)} \\rangle - (S_z + 1) (S_z + 2)}\n",
    "$$\n",
    "\n",
    "通过这种方式获得的 $\\langle S^2 \\rangle$ 近似值可以相当精确，比 $\\langle S^{2(0)} \\rangle + \\langle S^{2(1)} \\rangle$ 还要接近精确值 $3.75$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.749999998117527"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S2_0 + (S4SD - S2_0**2) / (S2_0 - (Sz + 1) * (Sz + 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "相信 $\\texttt{S4SD}$ 的存在意义是用于计算 Schlegel [^Schlegel-Schlegel.JCP.1986.84] 式 eq (24) 中的 $\\langle \\tilde \\Phi_1 | \\tilde \\Phi_1 \\rangle = \\langle \\Phi_0 | A_{s + 1}^\\dagger A_{s + 1} | \\Phi_0 \\rangle$；但关于这一关系我还不确定是否正确。后面计算 PMP2 能量时会使用上 $\\texttt{S4SD}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 自旋污染矫正能量计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `EPUHF` $E_\\mathrm{PUHF}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据 Schlegel [^Schlegel-Schlegel.JCP.1986.84] eq (22)，PUHF 能量可以表达为\n",
    "\n",
    "$$\n",
    "E_\\mathrm{PUHF} = E_\\mathrm{UHF} + \\frac{1}{\\langle \\Psi_0 | \\hat P_s | \\Psi_0 \\rangle} \\sum_{i \\bar j a \\bar b} \\langle \\Psi_0 | \\hat H | \\Psi_{i \\bar j}^{a \\bar b} \\rangle \\langle \\Psi_{i \\bar j}^{a \\bar b} | \\hat P_s | \\Psi_0 \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中，$\\hat P_s$ 算符称为 Löwdin 算符 [^Lowdin-Lowdin.PR.1955.97] eq (7)，\n",
    "\n",
    "$$\n",
    "\\hat P_s = \\prod_{k \\neq s}^{N / 2} \\frac{\\hat S^2 - k (k + 1)}{s (s + 1) - k (k + 1)}\n",
    "$$\n",
    "\n",
    "相当于将自旋不纯的波函数纯化为自旋量子数为 $s$ 的态。在实际使用中，通常使用 $\\hat A_{s + 1} \\simeq \\hat P_s$ 替代；关于这段讨论可以参考 Rossky and Karplus [^Rossky-Karplus.JCP.1980.73] section V.A 的讨论，而下面公式的形式参考 Schlegel [^Schlegel-Schlegel.JCP.1986.84] eq (14)；其中，$s$ 一般取 $S_z$：\n",
    "\n",
    "$$\n",
    "\\hat A_{s + 1} = \\frac{\\hat S^2 - (s + 1)(s + 2)}{\\langle S^{2(0)} \\rangle - (s + 1)(s + 2)}\n",
    "$$\n",
    "\n",
    "关于 $\\hat A_{s + 1}$，一个显然的性质是 $\\langle \\Psi_0 | \\hat A_{s + 1} | \\Psi_0 \\rangle = 1$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了程序方便，定义下述临时变量 `Y`\n",
    "\n",
    "$$\n",
    "Y = \\langle S^{2(0)} \\rangle - (S_z + 1) (S_z + 2)\n",
    "$$\n",
    "\n",
    "那么 `D_EPUHF`\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Delta E_\\mathrm{PUHF} &= \\sum_{i \\bar j a \\bar b} t_{i \\bar j}^{a \\bar b} D_{i \\bar j}^{a \\bar b} \\cdot \\langle \\Psi_{i \\bar j}^{a \\bar b} | \\frac{\\hat S^2}{Y} | \\Psi_0 \\rangle \\\\\n",
    "&= - \\frac{1}{Y} \\sum_{i \\bar j a \\bar b} t_{i \\bar j}^{a \\bar b} D_{i \\bar j}^{a \\bar b} S_{i \\bar b} S_{\\bar j a}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0010039336320381543"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = S2_0 - (Sz + 1) * (Sz + 2)\n",
    "D_EPUHF = - 1 / Y * (t2_ab * D2_ab * S_ia[:, None, None, :] * S_ai.T[None, :, :, None]).sum()\n",
    "D_EPUHF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因而 $E_\\mathrm{PUHF} = E_\\mathrm{UHF} + \\Delta E_\\mathrm{PUHF}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-73.04614628727792"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scf_eng.e_tot + D_EPUHF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian 的参考值是 -73.046146318。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `EPMP2` $E_\\mathrm{PMP2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据 Schlegel [^Schlegel-Schlegel.JCP.1986.84] eq (24)，PMP2 能量可以表达为\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Delta E_\\mathrm{PMP2} = \\Delta E_\\mathrm{PUHF} \\left( 1 - \\frac{\\langle \\Phi^{(1)} | \\hat A_{s + 1} | \\Psi_0 \\rangle}{\\langle \\Phi_0 | \\hat A_{s + 1}^\\dagger \\hat A_{s + 1} | \\Psi_0 \\rangle} \\right)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关于上式的分数项，分子部分可以写为\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle \\Phi^{(1)} | \\hat A_{s + 1} | \\Psi_0 \\rangle\n",
    "= \\langle \\Phi^{(1)} | \\frac{\\hat S^2}{Y} - \\frac{(s + 1)(s + 2)}{Y} | \\Psi_0 \\rangle = \\frac{1}{2} \\frac{\\langle S^{2(1)} \\rangle}{Y}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "而关于分子项，参考在 $\\texttt{S4SD}$ 的讨论，\n",
    "\n",
    "$$\n",
    "\\langle \\Phi_0 | \\hat A_{s + 1}^\\dagger \\hat A_{s + 1} | \\Psi_0 \\rangle \\simeq \\langle S^2 \\rangle - \\langle S^{2(0)} \\rangle = \\frac{\\texttt{S4SD} - \\langle S^{2(0)} \\rangle^2}{Y^2}\n",
    "$$\n",
    "\n",
    "但作者不能断定上述论断的正确性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将分子、分母的结果代入 $\\Delta E_\\mathrm{PMP2}$ 的算式中，可以得到 `D_EPMP2`\n",
    "\n",
    "$$\n",
    "\\Delta E_\\mathrm{PMP2} = \\Delta E_\\mathrm{PUHF} \\left( 1 - \\frac{1}{2} \\frac{\\langle S^{2(1)} \\rangle \\cdot Y}{\\texttt{S4SD} - \\langle S^{2(0)} \\rangle^2} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0005710125302116336"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_EPMP2 = D_EPUHF * (1 - 0.5 * S2_1 * Y / (S4SD - S2_0**2))\n",
    "D_EPMP2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因而 $E_\\mathrm{PMP2} = E_\\mathrm{UMP2} + \\Delta E_\\mathrm{PMP2}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-73.07218055803808"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp2_eng.e_tot + D_EPMP2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gaussian 的参考值是 -73.072180589。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "至此，我们已经完成了使用 PySCF 的功能与 NumPy 重复 Gaussian 的 PUHF、PMP2 的能量结果了。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 修订时间轴"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 2019/08/30 写完文档；文档基于 2019/08/13 的一份笔记。\n",
    "\n",
    "- 2019/09/01 补充一部分推导。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Chen-Schlegel.JCP.1994.101]: Chen, W.; Schlegel, H. B. Evaluation of S2 for Correlated Wave Functions and Spin Projection of Unrestricted Moller–Plesset Perturbation Theory. *J. Chem. Phys.* **1994**, *101* (7), 5957–5968. doi: [10.1063/1.467312](https://doi.org/10.1063/1.467312).\n",
    "\n",
    "[^Szabo-Ostlund.Dover.1996]: Szabo, A.; Ostlund, N. S. *Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry)*; Dover Publications, 1996.\n",
    "\n",
    "[^Schlegel-Schlegel.JCP.1986.84]: Schlegel, H. B. Potential Energy Curves Using Unrestricted Mo/ller–Plesset Perturbation Theory with Spin Annihilation. *J. Chem. Phys.* **1986**, *84* (8), 4530–4534. doi: [10.1063/1.450026](https://doi.org/10.1063/1.450026).\n",
    "\n",
    "[^Amos-Hall.PRSLA.1961.263]: Amos, A. T.; Hall, G. G. Single Determinant Wave Functions. *Proc. R. Soc. Lond. A* **1961**, *263* (1315), 483–493. doi: [10.1098/rspa.1961.0175](https://doi.org/10.1098/rspa.1961.0175).\n",
    "\n",
    "[^Lowdin-Lowdin.PR.1955.97]: Lowdin, P.-O. Quantum Theory of Many-Particle Systems. III. Extension of the Hartree-Fock Scheme to Include Degenerate Systems and Correlation Effects. *Phys. Rev.* **1955**, *97* (6), 1509–1520. [10.1103/physrev.97.1509](https://doi.org/10.1103/physrev.97.1509).\n",
    "\n",
    "[^Rossky-Karplus.JCP.1980.73]: Rossky, P. J.; Karplus, M. Spin Dependent Properties of Perturbed Wave Functions: An Analytic Comparison of the Exact, UHF, and Spin-Projected UHF States. *J. Chem. Phys.* **1980**, *73* (12), 6196–6214. [10.1063/1.440115](https://doi.org/10.1063/1.440115)."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "224px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
